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ABSTRACT 

Hyperion is a new three-dimensional dust continuum Monte-Carlo radiative transfer code that is designed to be as 
generic as possible, allowing radiative transfer to be computed through a variety of three-dimensional grids. The main 
part of the code is problem-independent, and only requires an arbitrary three-dimensional density structure, dust 
properties, the position and properties of the illuminating sources, and parameters controlling the running and output 
of the code. Hyperion is parallelized, and is shown to scale well to thousands of processes. Two common benchmark 
models for protoplanetary disks were computed, and the results are found to be in excellent agreement with those 
from other codes. Finally, to demonstrate the capabilities of the code, dust temperatures, SEDs, and synthetic multi- 
wavelength images were computed for a dynamical simulation of a low-mass star formation region. Hyperion is being 
actively developed to include new features, and is publicly available (http://www.hyperion-rt.org). 
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1. Introduction 

The investigation of astrophysical sources via multi- 
wavelength observations requires understanding the trans- 
fer of radiation through dust and gas in order to reliably de- 
rive geometrical, physical, and chemical properties. While 
a small subset of problems can be solved analytically, most 
require a numerical solution. 

A variety of techniques have been developed to this 
effect, and until recently, most relied on a direct nu- 
merical solution to the differential equation of radiative 
transfer. Examples include 2"^^ order finite differences 
(Steinacker et al. 2003), 5"^ order Runge-Kutta integra- 
tion (Steinacker et al. 2006), and variable Eddington ten- 
sors (Dullemond & Turolla 2000; DuUemond 2002). 

While these techniques can often provide results accu- 
rate to arbitrary precision, and have been very successful for 
one-dimensional problems, they become increasingly com- 
plex when applied to two- and three-dimensional problems. 
Within the last decade, the Monte- Carlo technique applied 
to radiative transfer has become an increasingly popular al- 
ternative that is well suited to arbitrary three-dimensional 
density distributions. Rather than solving the equation of 
radiative transfer directly, this method relics on random 
sampling of probability distribution functions to propagate 
photon packets through a grid of constant density and tem- 
perature cells. 

In its simplest implementation, Monte-Carlo radiative 
transfer considers single-frequency photon packets that arc 
emitted by sources, and are propagated through a grid of 
constant density cells. Each photon packet travels a cer- 
tain optical depth, before being either scattered, or ab- 
sorbed and immediately re-emitted (to conserve energy). 
These photon packets can then continue to travel and in- 



teract until they escape the grid. The optical depth, type 
of interaction, and frequency of re-emitted photon pack- 
ets, as well as scattcring/rc-cmission directions arc all sam- 
pled from probability distribution functions. Photon pack- 
ets escaping from the grid can be used to compute synthetic 
observations, such as spectral energy distributions (SEDs) 
and images. By repeating this process for large numbers of 
photon packets, the signal-to-noise of the synthetic obser- 
vations can be increased. 

In cases where local thermodynamic equilibrium (LTE) 
can be assumed, the energy absorbed in each cell can be 
used to compute the equilibrium temperature. Since the 
emissivity depends on the temperature, the radiative trans- 
fer has to be run iteratively to obtain a self-consistent so- 
lution for the temperature. 

Various methods have been developed to improve the 
performance of Monte-Carlo radiative transfer codes, which 
would otherwise be inefficient in many cases. For exam- 
ple, rather than simply using discrete absorptions in cells 
to compute the equilibrium temperature, Lucy (1999) sug- 
gested summing the optical depth to absorption of photon 
packet paths through each grid cell, requiring far fewer pho- 
tons to achieve a comparable signal-to-noise in the derived 
temperatures. Bjorkman & Wood (2001) developed an al- 
gorithm whereby the temperature of a cell is updated imme- 
diately each time a photon is absorbed, and the frequency 
of the re-emitted photon is sampled from a difference emis- 
sivity function that introduces a correction to account for 
the fact that previous photons were sampled from emissivi- 
ties for different temperatures. This algorithm is sometimes 
referred to simply as 'immediate re-emission', but this is a 
misnomer since other algorithms also immediately re-emit 
a photon following an absorption to conserve energy (e.g. 
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Lucy 1999); instead, what this algorithm introduces is the 
immediate temperature correction aspect. 

For synthetic observations, the peeling-off technique 
(Yusef-Zadeh et al. 1984) greatly improves the signal-to- 
noise of SEDs and images: every time a photon packet is 
emitted, scattered, or re-emitted after an absorption, the 
probability that the emitted or scattered photon would 
reach the observer is used to build up the observations. 
Thus each photon contributes multiple times to the SEDs 
or images. Another example is that of forced first scatter- 
ing (e.g. Mattila 1970; Wood & Reynolds 1999), which can 
be used in optically thin cases to force photon packets to 
interact with the dust (this requires weighting the photons 
accordingly to avoid any bias due to the forcing). 

In standard Monte-Carlo, the number of interactions a 
photon packet will undergo before traveling a certain dis- 
tance increases approximately with the square of the den- 
sity, so that in very optically thick cells, photon packets may 
become effectively trapped. This problem can be avoided 
by locally making use of the diffusion approximation to 
solve the radiative transfer in these regions. For example, 
Min et al. (2009) developed a modified random walk pro- 
cedure based on the diffusion approximation that can dra- 
matically reduce the number of steps required for a photon 
to escape from a grid cell of arbitrarily large optical depth. 

A large collection of dust continuum Monte-Carlo 
radiative transfer codes has been developed (e.g. 
Wolf et al. 1999; Harries 2000; Gordon et al. 2001; 
Misselt et al. 2001; Wood et al. 2001, 2002a,b; Wolf 
2003; Stamatellos & Whitworth 2003; Whitney et al. 
2003a,b, 2004; DuUemond & Dominik 2004; Jonsson 
2006; Pinte et al. 2006; Min et al. 2009), each including 
some or all of the above optimizations, as well as other 
optimizations not mentioned here. While some of the early 
codes assumed spherical or axis-symmetric geometries for 
simplicity, many have since been adapted to compute fully 
arbitrary three-dimensional distributions. In addition to 
dust continuum radiative transfer, some codes can also 
compute non-LTE line transfer (Carciofi & Bjorkman 
2006, 2008), or photoionization (e.g. Ercolano et al. 2003, 
2005, 2008). 

This paper presents Hyperion, a new dust-continuum 
Monte-Carlo radiative transfer code that is designed to be 
applicable to a wide range of problems. Hyperion imple- 
ments many of the recent optimizations to the Monte- Carlo 
technique discussed above, was written from the start to 
be a parallel code that can scale to thousands of processes, 
and is written in a modular and extensible way so as to be 
easily improved in future. It can treat the emission from 
an arbitrary number of sources, can include multiple dust 
types, and can compute the anisotropic scattering of polar- 
ized radiation using fully numerical scattering phase func- 
tions. It uses the Lucy (1999) iterative method to determine 
the radiative equilibrium temperature, but does not use 
the Bjorkman & Wood (2001) temperature correction tech- 
nique, as the latter is much more challenging to parallelize 
efficiently. Thanks to the modular nature of the code, the 
radiative transfer can be computed on a number of differ- 
ent three-dimensional grid types, and additional grid types 
can be added in future. Hyperion can compute SEDs and 
multi-wavelength images and polarization maps. The code 
is released under an open source license, and is hosted on 
a service that allows members of the community to easily 
contribute to the code and documentation. 



Section 2 gives an overview of the implementation of the 
code. Section 3 discusses the efhciency of the parallelized 
code. Section 4 presents the results for two benchmark mod- 
els of protoplanetary disks. Finally, Section 5 demonstrates 
the capabilities of the code by computing temperatures, 
SEDs, and synthetic images for a simulation of a star- 
formation region. The availability of the code and plans 
for the future are discussed in Section 6. 

2. Code overview 

The code is split into two main components. The first, 
which carries out the core of the radiative transfer calcu- 
lation, is implemented in Fortran 95/2003 for high perfor- 
mance. This part of the code is problem-independent: the 
input (bundled into a single file) consists of an arbitrary 
three-dimensional density structure as well as dust proper- 
ties, a list of sources, and output parameters. This input is 
used by the Monte-Carlo radiative transfer code to compute 
temperatures, SEDs, and images. Therefore, it is possible to 
use either gridded analytical density structures, or arbitrary 
density grids from simulations. At the moment, Hyperion 
supports several types of three-dimensional grids (§2.1.2) 
and the modular nature of the code will make it easy to add 
support for additional grid types in the future. It is possi- 
ble to specify an arbitrary number of dust types (within 
computational limits), which allows models to have differ- 
ent effective grain size distributions and compositions in 
different grid cells. 

The second component of the code consists of an object- 
oriented Python library that makes it easy to set up the 
input file for arbitrary problems from a single script. This 
library includes pre-defined analytical density structures for 
common problems such as flared disks and rotationally flat- 
tened envelopes and will also include scripts to import den- 
sity structures from simulations. Post-processing tools are 
also provided in the Python library to analyze the results 
of radiative transfer models. 

The present section describes the algorithm for the main 
radiative transfer code. The code first reads in the inputs 
(§2.1), then propagates photon packets through the grid 
(§2.2) for multiple iterations to compute the energy ab- 
sorbed in each cell (§2.3). Once the absorbed energy calcu- 
lation has converged (§2.4), the code computes SEDs and 
images (§2.6). 

2.1. Inputs to the code 
2.1.1. Sources 

A model can include any number of sources of emission 
(within computational limits). Each source is characterized 
by a bolometric luminosity, and the frequencies of the emit- 
ted photon packets are randomly sampled such that the 
emergent frequency distribution of the packets reproduces 
a user-defined spectrum. The total number of photons to 
emit from sources is set by the user. A number of differ- 
ent source types can be used - at the moment, the code 
supports: 

— Isotropic point sources. 

— Spherical sources with or without limb darkening. These 
can include arbitrary numbers of cool or hot spots, each 
with different positions, sizes, and with their own spec- 
trum. 



2 



Robitaille: Hyperion 



— Diffuse sources where flux is emitted from within grid 
cells according to a three-dimensional probability distri- 
bution function (this can be useful for unresolved stellar 
populations in galaxy models, or for accretion luminos- 
ity emitted via viscous energy dissipation in protoplan- 
etary disks). 

— External isotropic sources, which can be used to simu- 
late an interstellar or intergalactic radiation field. 

Each photon packet emitted is characterized by a position, 
direction vector, frequency, and a Stokes vector (I, Q, U, V) 
that describes the total intensity and the linear and circular 
polarization. 



2.1.2. Dust density grid 

The code is written in a modular way that allows sup- 
port for different grid geometries to be easily added. At the 
moment, three-dimensional cartesian, spherical-polar, and 
cylindrical-polar grids can be used, as well as two types of 
adaptive cartesian grids. The first is a standard octree grid, 
in which each cubic cell can be recursively split equally 
into eight smaller cubic cells. The second is the type of 
grid commonly used in adaptive mesh refinement (AMR) 
hydrodynamical codes. Here, a coarse grid is first defined 
on the zero-th level of refinement, and with increasing lev- 
els, increasingly finer grids can be used in areas where high 
resolution is needed. Because they concentrate the resolu- 
tion where it is needed, variable resolution grids such as 
octrees and AMR allow radiation transfer to be computed 
on density grids with effective resolutions that would be 
prohibitive with regular cartesian grids. 



2.1.3. Dust properties 

The dust properties required are the frequency-dependent 
mass extinction coefficient Xu ^J^d albedo lo^^ as well as the 
scattering properties of the dust. At this time, Hyperion 
supports anisotropic wavelength-dependent scattering of 
randomly oriented grains, using a 4-element Mueller ma- 
trix (Chandrasekhar 1960; Code & Whitney 1995): 
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Support for aligned non-spherical dust grains, which are 
described by a full 16-element matrix, will be implemented 
in future (e.g. Whitney & Wolff 2002). 

To keep the Fortran code as general as possible, 
the mean opacities and emissivities of the dust are pre- 
computed by the Python library as a function of the spe- 
cific energy absorption rate of the dust rather than the dust 
temperature (c.f. §2.3). For dust in LTE, the emissivities 
are given by ~ k^, (T) , and the mean opacities are 
the usual Planck and Rosseland mean opacities to extinc- 
tion and absorption. However, it is also possible to specify 
mean opacities and emissivities that do not assume LTE 
(e.g. §2.5). Thus, assumptions about LTE are made at the 
level of the dust files, rather than in the radiative transfer 
code itself. 



2.2. Photon packet propagation 

The code implements the propagation of photon packets 
in the following way: a photon packet is emitted from a 
source, randomly selected from a probability distribution 
function defined by the relative luminosity of the different 
sources. This sampling can be done either in the standard 
way to give a number of photon packets proportional to 
the source luminosity, or to give equal numbers of pho- 
tons to each source, which requires weighting the energy 
of the photons. The direction and frequency of the photon 
packet are randomly sampled accordingly for the type and 
the spectrum of the source it originates from, using stan- 
dard sampling with no weighting. A random optical depth 
to extinction t is sampled from the probability distribution 
function exp (— r) by sampling a random number f uni- 
formly between and 1, and taking t = — In The photon 
packet is then propagated along a ray until it either escapes 
the grid, or reaches the sampled optical depth, whichever 
happens first. If the photon packet has not escaped the grid, 
it will then interact with the dust. A random number ^ is 
sampled uniformly between and 1, and if it is larger than 
the albedo of the dust, the photon packet is absorbed; oth- 
erwise it is scattered. Once the photon packet is scattered 
or rc-emittcd, a new optical depth is sampled, and the pro- 
cess is repeated until the photon packet escapes from the 
grid. 

Very optically thick regions are an issue in basic 
Monte-Carlo radiative transfer, as photon packets can get 
trapped in these regions and require in some cases mil- 
lions of absorptions/re-emissions and scatterings to escape. 
Recently, Min et al. (2009) proposed a modified random 
walk (MRW) algorithm that allows photon packets to prop- 
agate efficiently in very optically thick regions by grouping 
many scatterings and absorptions/re-emissions into single 
larger steps. The photon packet propagation described pre- 
viously is done in a grid made up of cells of constant den- 
sity and temperature. Therefore if the mean optical depth 
to the edge of a cell is much larger than unity, one can set 
up a sphere whose radius is smaller than the distance to 
the closest wall, inside which the density will be constant, 
and travel to the edge of a sphere in a single step using the 
diffusion approximation, thus avoiding the need to com- 
pute millions of interactions. Hyperion includes the im- 
plementation of the MRW algorithm described in Robitaille 
(2010). 

2.3. Temperature/Energy absorption rate calculation 

Hyperion uses the iterative continuous absorption method 
proposed by Lucy (1999): in the first iteration, the specific 
energy absorption rate of the dust is computed in each cell 
using 



where At is the time over which photon packets are emitted 
(taken to be Is in the code), V is the cell volume, e is the 
energy of a photon packet, i is the path length traveled - 
which depends on the density, as does the number of path 
lengths being added - and Ki, is the mass absorption coeffi- 
cient. The temperature T of the dust can be found from A 
by balancing absorbed and emitted energy, assuming LTE: 



ATrKp{T)B{T) ^ A 



(3) 
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where Kp{T) is the Planck mean mass absorption coeffi- 
cient, and B(T) = (tT/7r)T^ is the integral of the Planck 
function. As mentioned in §2.1.3, Hyperion docs not com- 
pute temperatures, but instead the mean opacities and 
emissivities of the dust, which are pre-computed, are tab- 
ulated as a function of A rather than T. In cases where 
the temperature is needed (for example if requested as out- 
put from the user), the temperature is computed on the fly 
using the pre-computed Planck mean opacities. 

For high optical depth problems, such as protoplanetary 
disks, some regions in the grid may see few or no photon 
packets, and reliable values for A or T can therefore not be 
directly computed. In this case, one can formally solve the 
diffusion approximation for the cells that see fewer than a 
given threshold of photon packets, using cells that do have 
reliable values as boundary conditions (see Appendix A). 
This is referred to as the partial diffusion approximation 
(PDA) in Min et al. (2009). 

2.4. Convergence 

The emissivity of the dust depends on A, which depends on 
the propagation of the photon packets and the frequency 
of the photon packets, which in turn depends on the emis- 
sivity of the dust, so one needs to compute the radiative 
transfer for several iterations before the values of A in each 
cell converge. A simple algorithm to determine convergence 
is included in the code. The main function used in the con- 
vergence algorithm is 

6 {xi ,X2) = max ( — , — ) (4) 
\X2 xi J 

This measures how different X2 is from xi and does not 
depend on the direction of the change. For example, a value 
of i5 = 2 means that the quantity has changed by a factor 
of 2. Because the change is expressed as a ratio, this means 
that large changes can be more easily expressed than using 
a simple fractional difference. 

At each iteration i, the code determines by how much 
the energy absorbed A has changed by computing for each 
cell j the change in the energy absorbed in each cell j, 

i?}^5(i}-\i}) (5) 

The quantile value at the p-th percentile of the Rj val- 
ues is then computed and compared to the value found 
during the previous iteration, Q^~^ . Finally, the change in 
this quantile is calculated using 

A' = S{Q'-\Q') (6) 

The specific energy absorption rate is considered to be con- 
verged once and A' have fallen below user specified val- 
ues Qthros and Athrcs- As an example, setting p = 99.9%, 
Qthres = 2., and Athrcs = 1-1 mcaus that convergence is 
achieved once 99.9% of the differences in specific energy 
absorption rates between iterations are less than a factor 
of 2, and once the 99.9% percentile value of the difference 
changes by less than a factor of 1.1 10%). 

Using this techniques is more robust to noise in the spe- 
cific energy values than simply requiring that all cells vary 
less than a given threshold, since it allows the user to ig- 
nore outliers by setting the percentile value appropriately. 



Of course, there is no guarantee that the criteria set by the 
user will be met at any iteration if the noise is too large, 
but in this case the user will see that there is not con- 
vergence, and can increase the number of photon packets, 
the maximum number of iterations, or use less stringent re- 
quirements for convergence. In any case, this method allows 
users to set quantitative requirements that are necessary to 
meet their scientific problem and that are more flexible than 
a single threshold on all values. 

2.5. PAH/VSG emission 

In addition to dust continuum radiative transfer, the code 
can also take into account any population where the opac- 
ity is independent of temperature or density, and for which 
the emissivity depends only on the specific energy absorbed 
(A) inside each cell (§2.1.3). While this does not permit ar- 
bitrary non-LTE radiative transfer, it does allow one to 
use an approximation of the emission from stochastically 
heated polycyclic aromatic hydrocarbons (PAHs) and very 
small non-thermal grains (VSGs) using a similar prescrip- 
tion to that given in Wood et al. (2008). 

The idea behind the algorithm is to pre-compute the 
average emissivity of an ensemble of PAHs and VSGs, as a 
function of the energy deposited by radiation into the PAHs 
and VSGs - the specific energy absorption rate A for the 
PAH and VSG populations - for a given irradiating spec- 
trum, and to then use these emissivities as look-up tables 
in the radiative transfer code. 

Since the opacities of PAHs and VSGs are typically 
strongly peaked in the UV and optical, this means that 
the UV and optical will dominate the energy exciting and 
being reprocessed by the PAHs and VSGs, and the emis- 
sivities are then being chosen based on the total strength 
of the UV and optical emission relative to the original tem- 
plate spectrum. The shape of the template spectrum does 
have a small impact on the pre-computed emissivities, but 
to first order, simply using the ratio of the total energy ab- 
sorbed should be adequate for estimating the importance 
of PAHs and VSGs as absorbers and emitters, and allows 
the impact on the SEDs and images to be studied. While 
the strength of the PAH features and VSG continuum in 
SEDs should be accurate to first order, the shape of the 
PAH features should be treated with caution. 

Note that the implementation discussed here differs in 
one respect from the Wood et al. (2008) prescription: the 
latter uses the mean intensity, rather than the energy inter- 
cepted by the PAHs and VSGs, to determine which emis- 
sivity file to use. Using the energy absorbed by the PAHs 
and VSGs may be more appropriate than using the mean 
intensity, since the latter would predict the same excitation 
for a fixed mean intensity, whether this intensity peaked in 
the UV or the millimeter, while using A would give a higher 
excitation for a spectrum peaking in the UV compared to 
one peaking in the mm. 

2.6. SEDs and images 

Once the specific energy absorption rate calculation has 
converged, SEDs and images can be computed. There are 
several methods to do this, and all of the following are 
implemented in the code. 
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2.6.1. Photon binning 

The easiest and most inefficient method to compute SEDs 
and images is to propagate the photon packets as for the 
initial iterations, and bin them all into viewing angles as 
they escape from the grid. This is very inefhcient, because 
each photon packet only contributes once to the SEDs and 
images, and only to one viewing angle. Furthermore, the 
viewing angle bins cannot be arbitrarily small, and there- 
fore the SEDs and images resulting from this are averaged 
over viewing angle. 

2.6.2. Peeling-off 

Yusef-Zadeh et al. (1984) introduced the concept of 
peeling-off, whereby at each scattering or re-emission, the 
probability p of the photon packet being scattered or re- 
emitted towards the observer immediately after the in- 
teraction is computed, and a photon packet with weight 
p exp (— r) is added to the SED and images, where r is 
the optical depth to reach the observer from the interac- 
tion. This results in much higher signal-to-noise than pho- 
ton packet binning, because each photon packet contributes 
several times to the SEDs and images at each viewing angle. 

2.6.3. Raytracing 

By far the most efficient method of computing SEDs and 
images is raytracing, which essentially consists of deter- 
mining the source function at each position in the grid, 
and solving - in post-processing - the equation of radia- 
tive transfer along lines of sight to the observer through 
the dust geometry. For thermal emission, this is relatively 
straightforward, because the source function in each cell 
is simply related to the mass and temperature or energy 
in the cell. For scattered light, unless the scattering phase 
function is isotropic, one needs to retain information about 
the angular and frequency dependence of the incident or 
scattered light. As discussed in Pinte et al. (2009), this is 
either computationally very expensive in terms of memory, 
or results in a loss of accuracy for strongly peaked scat- 
tering phase functions. In the current implementation of 
Hyperion, raytracing can be used for the source and dust 
emission, allowing excellent signal-to-noise to be achieved 
at long wavelengths where traditional Monte-Carlo only 
produces very few photon packets. Raytracing for scattered 
light will be implemented as an option in future since ~ if 
adequate computational resources are available - it can pro- 
vide excellent signal-to-noise significantly faster than con- 
ventional peeling-off as for source and dust emission. 

2.6.4. Monochromatic radiative transfer 

The default mode of the binning (§2.6.1), peeling-off 
(§2.6.2), and raytracing (§2.6.3) algorithms is to produce 
SEDs and images with finite width wavelength/frequency 
bins. However, in some cases it is desirable to compute 
SEDs or images at exact wavelengths/frequencies. When 
this is the case, the peeling-off algorithm for computing the 
scattered light contribution to the SED or images has to 
be modified. In this case, the scattered light is separated 
into two contributions, namely the scattered light from the 
sources, and the scattered dust emission. To compute the 
scattered light from the sources, photon packets are emit- 



ted by all the sources at the fixed wavelengths / frequencies 
required. Each time a photon packet scatters, a photon 
packet is peeled-off, and each time a photon packet is ab- 
sorbed, it is terminated. Thus, there is no immediate re- 
emission following an absorption. The next step is to com- 
pute the scattered dust emission. The specific energy ab- 
sorption rate is used to compute the emissivity at the fixed 
wavelengths/frequencies inside each cell. To emit a photon 
packet, a random cell is selected in the grid, and the pho- 
ton packet is emitted randomly within the cell, carrying an 
amount of energy proportional to the local emissivity. The 
photon packet is then propagated, and a photon packet is 
peeled-off at each scattering. As before, the photon packet 
is terminated once it is absorbed. This algorithm ensures 
conservation of the total scattered light contribution. 

2. 7. Additional user options 

2.7.1. Uncertainties 

When computing SEDs and images, Hyperion allows the 
user to compute uncertainties, which uses the scatter in the 
photon packet flux values to derive errors in the flux at each 
wavelength/frequency and in each aperture or pixel. The 
fact that a given SED or image can be constructed using a 
combination of techniques, such as peeling-off and raytrac- 
ing - which can produce very different signal-to-noise - is 
properly taken into account by computing the uncertainties 
for the contribution from each technique to the final SEDs 
or images separately and then combining them. 

2.7.2. Photon tracking 

Hyperion offers the option for the user to track the origin 
of each photon packet to split SEDs and images into dif- 
ferent components. A basic mode allows SEDs and images 
to be split into contributions from sources and from dust, 
while a more detailed mode allows the flux to be split into 
individual sources and dust types. In both cases, the flux 
can be split further into photon packets reaching the ob- 
server directly, and photon packets having been scattered 
since the last emission/re-emission and before reaching the 
observer. 

2.7.3. Dust sublimation 

Users have the option to specify a dust sublimation spe- 
cific energy absorption rate for each dust type, and three 
different dust sublimation modes are possible: 

— the specific energy absorption rate can simply be capped 
to the maximum value, without changing the density. 

— the dust can be completely removed from cells exceeding 
the maximum specific energy absorption rate. 

— the dust density can be reduced but not set to zero. 
This can be useful because in optically thick cells, if ra- 
diation originates from the cell (for example luminosity 
from viscous dissipation) , the specific energy absorption 
rate can be high because the radiation is trapped. If the 
density had been lower but non-zero, the specific energy 
absorption rate might not have exceeded the maximum 
specified. In this situation, the dust density should be 
reduced but not set to zero. 
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Fig. 1. The main panels show the wall clock speedup relative to serial execution for a model of a protoplanetary disk. The 
right panel shows the same model as that in the left panel, but with 10 times more photon packets for every iteration. 
The open circles show the average values of 10 executions of the code, and the uncertainties in these values derived from 
the scatter are not shown as they are smaller than the data points. The dashed curves show the speedup expected from 
a perfectly parallelized code. The solid curves show the best fit to the speedup curve assuming Amdahl's law (described 
in text). The gray shaded areas around the curves show the 68.3% confidence interval. The maximum speedups derived 
from the best fits are shown as horizontal dotted lines. The inset panels show the probability of different P values, with 
the 68.3% confidence intervals shown in shaded gray. 



2.7.4. Forced first scattering 

As mentioned in Section 1, forced first scattering (e.g. 
Mattila 1970; Wood & Reynolds 1999) can be used to im- 
prove the signal-to-noisc of scattered radiation in optically 
thin dust. HYPERION includes an implementation of this 
algorithm. 



3. Parallelized performance 

The reason for not using the Bjorkman & Wood (2001) 
immediate temperature correction method in Hyperion 
is that each time a photon gets absorbed and re-emitted, 
the specific energy absorption rate grid has to be updated. 
Thus, it is not possible to compute the propagation of mul- 
tiple photons simultaneously and to then combine the re- 
sults, and codes using this method can therefore not eas- 
ily be parallelized. By simply using an iterative approach 
to computing the specific energy absorption rate using the 
Lucy (1999) continuous absorption method, one has the 
advantage that within an iteration, the problem is "em- 
barrassingly parallel" . Each process can propagate photon 
packets independently, and at the end of the iteration, the 
energy absorbed in each cell is synchronized over processes. 
Similarly, when SEDs, images, or polarization maps are 
computed, these can be computed by separate processes, 
and combined at the end of the calculation. 

The code has been parallelized using the Message 
Passing Interface (MPI) . Since different cores or nodes can 
have inhomogeneous performance, rather than dividing the 
total number of photon packets equally between processes, 
the task is divided into chunks of photon packets that arc 
small enough that there are many more chunks than to- 
tal number of processes, but large enough that the pro- 
cesses do not communicate too often. Each process then 
computes one batch of photon packets, and reports back 
to the main process to find out whether to process an- 



other batch or whether to stop. This incurs a very small 
overhead, since the request consists essentially of a single 
number both ways. Only at the end of the iteration, once 
all processes have received the signal to stop, are the re- 
sults combined. This last step incurs an overhead, which 
scales with the number of processes and the grid, image, 
or SED resolution. However, in most cases, the overhead is 
negligible compared to the main computation. 

The parallel efficiency of the code depends on the par- 
ticular model being computed and the number of photon 
packets requested. One way to quantify the speedup for a 
given model is to compare the fraction of time spent in the 
serial execution of the code, such as the startup phase where 
the data is read in from the input file, or the time between 
iterations, to the fraction of time spent in the "embarrass- 
ingly parallel" photon packet propagation. If one writes the 
fraction of the code that is truly parallel as P, then the 
speedup obtained by running the code on N parallel pro- 
cesses is: 



speedup = 



1 



(1 - P) + P/N 



(7) 



This is often referred to as Amdahl's law (Amdahl 1967). 
One consequence of Equation (7) is that as TV tends to 
infinity, the speedup tends to a finite maximum 



speedup,,,^^ = 



1 



(8) 



In this case, the parallelized part of the code runs infinitely 
fast, and the execution time is given purely by the serial 
portion of the code. 

The value of P is dependent on the choice of the model 
being computed and number of photon packets. This is 
demonstrated here with a model of a protoplanetary disk 
(taken from §4.2) computed on different numbers of cores, 
ranging from iV = 1 to = 512 in powers of 2. The wall 
clock times of the parallel computations are compared to 
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Radius (AU) 6 (degrees) 

Fig. 2. Temperature results for the Paseucci ct al. (2004) disk benchmark, for the Ty = 100 model. The left panel shows 
a radial temperature profile close to the mid-plane {9 ~ 2.5°), both in absolute terms (top) and relative to the reference 
code RADICAL (bottom) . The right panel shows a vertical cut through the disk at a cylindrical radius of 2 AU. The 
black lines show the results from Hyperion, while the gray lines show the results from the codes tested in Paseucci et al. 
(2004) 



the true serial runtime rather than the parallel version run- 
ning with one process, since these will have slightly differ- 
ent runtimes. The times used are wall clock times^, rather 
than CPU times. For each number of cores, the model was 
run ten times to obtain a mean speedup and a measure 
of the scatter from one run to the next. The results are 
shown in Figure 1 for two instances of the model - one 
with ten times more photon packets than the other. The 
models were run on the Harvard Odyssey cluster, allow- 
ing N up to 512. For the model with the lower number 
of photon packets^, the speedups obtained are well fit by 
P = 0.99658 ± 0.00014, while the speedups for the model 
with the higher number of photon packets were well fit by 
P = 0.999373 ± 0.00006 where the uncertainties were de- 
rived as a formal confidence interval using the surface. 
The reduced values for the fits were 1.01 and 3.25 re- 
spectively, indicating good fits. The P values translate into 
maximum speedups of 292.6l}?;| and 1604:At\l^i respec- 
tively. These values differ by almost an order of magnitude, 
which is expected, since the only difference between the two 
models is that the CPU time of the parallel section of the 
code was increased by a factor of ten, while the serial part 
of the code remained identical (same input / output and syn- 
chronization overheads). In other words, both models tend 
to the same runtime, since the serial portion of the code is 
the same in both cases. 

The speedup values presented here are purely illustra- 
tive, and should not be assumed to hold for all models, but 
they serve to show that the efficiency of the code does have 
a finite limit which depends on the model set-up. As for 



^ This can also be referred to as 'real' time 

^ For the model with the lower number of photon packets, the 
results for A*' = 512 were not used, since the time to start up 
the processes on the nodes dominated the wall clock time, and 
therefore the speedup values obtained were a measure of the 
efficiency of the cluster than of the code 



most parallel codes, beyond a certain number of parallel 
processes, the wall clock runtime is completely dominated 
by the time to execute the serial part of the code (typically 
input/output and calculations between iterations such as 
the PDA). Users are encouraged to determine the optimal 
number of processes for the problem they are considering. 



4. Benchmarks 

4.1. Paseucci et al. (2004) benchmark 
4.1.1. Definition 

Paseucci et al. (2004) defined a benchmark problem for 2D 
radiative transfer, consisting of a flared disk around a point 
source. The central source is a point source with a 5,800 K 
blackbody spectrum, and a luminosity of 1.016 Lq . The disk 
extends from 1 to 1,000 AU, and has a density structure 
given by: 



p{r,z) = Pa 



ro 



exp 



ho {r/roY 



(9) 



where r is the cylindrical polar radius, rg = 500 AU, 
ho = 125 • \/27^AU, a = 1., and /3 = 1.125. The disk is 
made up of dust grains with a single size of 0.12 fim, a den- 
sity of 3.6g/cm^, and composed of astronomical silicates 
with properties given by Draine & Lee (1984). Scattering 
is assumed to be isotropic. 

The benchmark consists of four cases, with visual (A = 
550 nm) optical depths, measured radially along the disk 
mid-plane, of ry =0.1, 1, 10, and 100. These correspond 
to disk (dust) masses of 1.1138 x lO"'^, 1.1138 x 10"^, 
1.1138 X 10-^ and 1.1138 x IO-^'Mq respectively. The cor- 
responding maximum optical depths perpendicular to the 
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Fig. 3. SED results for the Pascucci et al. (2004) disk benchmark, for the four disk masses used. In eaeh panel, the top 
section shows the SED obtained using Hyperion (black line) , compared to the reference code from Pascucci et al. (2004) 
(RADICAL; gray circles), for the three viewing angles used (12.5, 42.5, and 77.5°). The dashed line shows the blackbody 
spectrum of the central source. The bottom three sections in each panel show the fractional difference between Hyperion 
and RADICAL (black line), and between the other codes in Pascucci et al. (2004) and RADICAL (gray lines). 



disk are of the order of ry =0.01, 0.1, 1, and 10 respec- 
tively. SEDs are computed at three viewing angles in each 
case: 12.5, 42.5, and 77.5°. 



4.1.2. Results 

The benchmark models were computed using a spheri- 
cal polar grid with dimensions (n^, ng, fi^) = (499,399,1) 
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Fig. 4. Temperature results for the Pinte et al. (2009) disk benchmark. The top-left and top-right panels shows the radial 
temperature profile at the disk mid-plane for the lowest and highest disk mass models. The bottom panels shows two 
vertical cuts through the disk for the highest disk mass model, at cylindrical radii of 0.2 and 200 AU respectively. In 
each panel, the top section shows the temperature profile in absolute terms, while the bottom section shows the results 
relative to the reference result from Pinte et al. (2009), which is the average of the results from MCFOST, MCMAX, 
and TORUS. The black lines show the results from Hyperion, while the gray lines show the results from MCFOST, 
MCMAX, and TORUS. 



and extending out to 1,300 AU. The SEDs were com- 
puted at specific wavelengths using the monochromatic 
radiative transfer described in §2.6. The MieX code 
(Wolf & Voshchinnikov 2004) was used to compute the op- 
tical properties of the dust. The convergence algorithm dis- 
cussed in §2.4 was used, with p = 99%, Qthros = 2, and 
Athres = 1.02. The models were run with enough photon 
packets to provide very high signal-to-noise in the temper- 
atures and SEDs'^, and were run using the parallel version 
of the code on a single machine with 8 cores, requiring a 
wall clock time of under 20 minutes for each model. The 
TV =0.1 and 1 models converged after 3 iterations, while 
the Ty =10 and 100 models converged after 4 iterations. 

Figure 2 shows for the most optically thick case (ry = 
100) the temperature profile for a fixed polar angle {6 = 
2.5°) as a function of radius r, and the temperature profile 



^ 10^ photons for each specific energy iteration, 10''' photons 
per wavelength for the monochromatic peeling-ofT for scattering 
for both the source and dust emission, and 10^ photons for the 
source and dust emission raytracing. 



for a fixed radius (r = 2 AU) as a function of polar angle 
9. The temperatures found by Hyperion are within the 
dispersion of the results from the other codes. 

Figures 3 shows the SEDs for the Hyperion code and 
the reference code in Pascucci et al. (RADICAL) for the 
four disk masses and three viewing angles, and the frac- 
tional difference between the SEDs is shown in each case. 
Also shown in gray are the differences for other codes pre- 
sented in Pascucci et al.. The SEDs from Hyperion are 
within the dispersion of results from the other codes. 

4.2. Pinte et al. (2009) benchmark 
4.2.1. Definition 

While the Pascucci et al. (2004) benchmark includes a case 
with fairly high visual optical depth through the mid-plane 
of the disk, the optical depths through realistic protoplan- 
etary disks are much higher (ry > 10^ in some cases), and 
are such that computing temperatures in the disk mid-plane 
becomes computationally challenging and requires approxi- 
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Fig. 5. SED results for the Pinte et al. (2009) disk benchmark, for the four disk masses used. In each panel, the top 
section shows the SED obtained using Hyperion (black line), compared to the average result of the codes used in 
Pinte et al. (2009) (gray circles), for four of the viewing angles (18.2, 75.5, 81.4, 87.1°). The dashed line shows the 
blackbody spectrum of the central source. The bottom four sections in each panel show the fractional difference between 
Hyperion and the reference result (black line), and between the other codes compared in Pinte et al. (2009) and the 
reference result (gray lines). 



mations to be made (e.g. §2.2 and §2.3). Pinte et al. (2009) of more massive disks with a more extreme radial density 
developed a complementary benchmark problem consisting profile to test codes in the limit of very optically thick disks. 
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Pixel offset 

Fig. 6. Image results for the Pinte et al. (2009) disk benchmark. The top two rows show the results for the i = 69.5° 
model, while the bottom two rows show the results for the i ~ 87.1° model. For each viewing angle, the resulting image 
is shown on the left on a power-law stretch (with power 1/4) , while the remaining plots show various cuts (indicated on 
the images), in absolute flux units, as well as relative to the reference result in Pinte et al. (2009), which is the average of 
the MCFOST, MCMAX, TORUS, and Pinball codes. The black hues show the results from Hyperion, while the gray 
lines show the resuhs from MCFOST, MCMAX, TORUS, and Pinball. As in Pinte et al. (2009), the cuts are 11 pixels 
wide to improve the signal-to-noise of the profiles. 



In addition, the benchmark also tests the ability to com- 
pute anisotropic scattering, and to reproduce images and 
polarization maps at various viewing angles. 

The benchmark case consists once again of a central 
star surrounded by a disk. The central star has a radius of 
2Rq and a temperature of 4000 K, with a blackbody spec- 
trum. The disk extends from 0.1 to 400 AU (with cylindri- 
cal edges), and the density is given by Equation (9), with 



ro = 100 AU, ho = 10 AU, a = 2.625, and f3 = 1.125. The 
disk is made up of dust grains with a single size of 1 /zm, a 
density of 3 . 5 g/cm^ , and composed of astronomical silicates 
with properties given by Weingartner & Draine (2001). 

The benchmark consists of four cases, with disk (dust) 
masses of 3 x 10'^ 3 x 10"'^, 3 x 10"^ and 3 x lO-^M©. 
While the disk masses are comparable to that of the 
Pascucci et al. benchmark problem, the smaller inner ra- 
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Fig. 7. Linear polarization map results for the Pinte et al. (2009) disk benchmark. The top two rows show the results 
for the i = 69.5° model, while the bottom two rows show the results for the i = 87.1° model. For each viewing angle, 
the resulting linear polarization map is shown on the left on a linear stretch, while the remaining plots show various cuts 
(indicated on the maps, and identical to those in Figure 6), both in absolute terms, and relative to the reference result in 
Pinte et al. (2009), which is the average of the MCFOST, MCMAX, and PinbaU codes. The black lines show the results 
from Hyperion, while the gray lines show the results from MCFOST, MCMAX, and PinbaU. As in Pinte et al. (2009), 
the cuts are 11 pixels wide to improve the signal-to-noise of the profiles. 



dius (0.1 vs 1 AU) and the much steeper surface density 
profile make this a more computationally challenging prob- 
lem. In fact, the I-band (A = 810 nm) optical depths, mea- 
sured radially along the disk mid-plane, are of the order of 
T/ = 10'^, lO'', 10^, and 10^ for the four disk masses respec- 
tively, orders of magnitude higher than for the Pascucci et 
al. benchmark. The corresponding maximum optical depths 



perpendicular to the disk arc of the order of t/ = 100, 10'^, 
lO"', and 10^ respectively. 

The benchmark test compares SEDs, images and polar- 
ization maps at 10 viewing angles corresponding to cos i = 
0.05 to 0.95 in steps of 0.10, corresponding to viewing angles 
from 18.2 (almost pole-on) to 87.1° (almost edge-on). The 
images and polarization maps are computed at 1 /im, and 
span 251 x 251 pixels, and a physical size of 900 x 900 AU. 
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Fig. 8. Three-color unconvolved synthetic images of the simulation, from the near-infrared (top left) and mid-infrared (top 
right) to the far- infrared (bottom left and right). The wavelengths and colors used are 1.25 /xm (blue), 1.65 //m (green), 
and 2.15 /xm (red) for the JHK image, 3.6 /im (blue), 4.5 /im (green), and 8.0 /im (red) for the Spitzer/IRAC image, 
70 /xm (blue), 110/xm (green), and 170 /im (red) for the Herschel/PACS image, and 250 /xm (blue), 350 /xm (green), 
and 500 /xm (red) for the Herschel /SPIRE image. All images are shown on an arcsinh stretch which allows a larger 
dynamic range to be shown by reducing the intensity of the brightest regions. The exact transformation is given by 
Vmax ■ asinh (mw/timax) /asinh(m) where v is the original flux, m = 30 is a parameter controlling the compression 
of the dynamic range, and Umax is 6, 4, and 2MJy/sr for J, H, and K; 2, 2, and 4MJy/sr for 3.6/xm, 4.5/xm, and 
8.0 /xm; 3000MJy/sr for the three PACS bands; and 1000, 500, and 250MJy/sr for SPIRE 250 /xm, 350 /xm, and 500 fM 
respectively. 
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Fig. 9. Three-color synthetic images of the simulation, from the near-infrared (top left) and mid-infrared (top right) 
to the far-infrared (bottom left and right). The images are binned to the pixel resolutions appropriate for ground 
based observations (JHK) and for Spitzer and Herschel observations respectively. The images were convolved with the 
appropriate instrumental PSFs, and include Gaussian noise with realistic levels. The wavelengths, stretches, and levels 
used for the images arc identical to those used in Figure 8. 
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Due to the single grain size, the scattering phase function 
oscillates strongly around 1 //m and at shorter wavelengths. 
Thus, these oscillations have a direct impact on scattered- 
light images computed at these wavelengths. While a single 
grain size, and therefore these oscillations, would not be 
seen in a real disk, this allows codes to test whether they 
correctly compute the scattering of photon packets. 

4.2.2. Results 

The benchmark models were computed using a cylindri- 
cal polar grid with dimensions {nr,nz,n^) = (499,399,1) 
extending out to the outer disk radius. The SEDs were 
computed at exact wavelengths using the monochromatic 
radiative transfer described in §2.6. As before, the MieX 
code was used to compute the optical properties of the dust. 
The convergence algorithm discussed in §2.4 was used, with 
p = 99%, Qthrcs = 2, and Athros = 1-02. The MRW (with 
7 = 2) and PDA approximations were used due to the high 
optical depths. The models were run with enough photon 
packets to provide very high signal-to-noise in the temper- 
atures, SEDs, and images'*^, and were run using the parallel 
version of the code on 32 cores, requiring a wall clock time 
of 12 to 17 hours to compute the temperature and SEDs 
for the Tj = 10'^ to 10^ models, and an additional 11 hours 
for the images and polarization maps. The temperatures 
converged after 7, 9, 9, and 10 iterations for the tj = 10'^, 
10*, 10^, and 10® models respectively. 

Figure 4 compares the mid-plane temperature profile in 
two of the models (r/ = 10"^ and t/ = 10®), as well as 
two vertical cuts in the most optically thick model (t/ = 
10®). Figure 5 compares the SEDs for the four models at 
four of the viewing angles (18.2, 75.5, 81.4, 87.1°). Finally, 
Figures 6 and 7 compare the total intensity and polarization 
fraction maps for the most optically thick model {tj = 10®). 
In all cases, the results produced by Hyperion are within 
the dispersion of results obtained from the other codes. 

5. Case study: simulated observations of a 
low-mass star-forming region 

5.1. The simulation 

To demonstrate the capabilities of the code, a radiative 
transfer calculation on a simulation of a low-mass star form- 
ing cloud was carried out in order to produce synthetic ob- 
servations from near-infrared to sub-milimeter wavelengths. 

The demonstration uses a simulation of low-mass 
star formation computed with the ORION AMR 
three-dimensional gravito-radiation-hydrodyamics code 
(Offner et al. 2009, and references therein). The simulation 
contains 185 of gas and dust in a box of side 0.65 pc, 
and resolves size-scales down to 32 AU, with an effective 
resolution of 4096'^. Offner et al. studied the effects of 
heating from protostars on the formation of stars, and 
therefore ran a simulation including radiative heating, and 
a control case with no radiation heating. The simulation 
used here is a snapshot of the one including radiative 
heating, for t iff , where iff is the free-fall time. 

* 10* photons for each specific energy iteration, 10*^ and 2 x 10*^ 
photons per wavelength for the monochromatic peeling-oflt for 
scattering of source and dust emission respectively (for SEDs; 
100 times more for the images and polarization maps), and 10^ 
photons for the raytracing for both the source and dust emission. 



Table 1. Noise levels added to the synthetic images in 
Figure 9. 



Band 


Noise a 




(MJy/sr) 


,J 


0.03 


H 


0.03 


K 


0.03 


IRAC 3.6 urn 


0.005 


IRAC 4.5 


0.005 


IRAC 8.0 


0.03 


PACS 70 fim 


20. 


PACS 110 //m 


30. 


PACS 170 /jm 


40. 


SPIRE 250 urn 


15. 


SPIRE 350 ^im 


10. 


SPIRE 500 nm 


5. 



5.2. Tlie radiative transfer model 

The radiative transfer through the star-forming region was 
simulated by taking into account both the luminosity from 
the forming stars (c.f. Appendix B of Offner et al. 2009), 
represented by sink particles in the simulation, and the 
interstellar radiation field at the solar neighborhood from 
Porter & Strong (2005), which includes contributions from 
the stellar emission, PAH emission, and far-infrared ther- 
mal emission^. The simulation was embedded at the cen- 
ter of a cube of side 1.95 pc and a uniform density of 
lO"^'^ g/cm to simulate an ambient medium. The dust 
properties were taken from Draine (2003a,b) using the 
Weingartner & Draine (2001) Milky Way grain size distri- 
bution A for i?v=5.5 and C/H = 30ppm renormalized to 
C/H = 42.6 ppm. In order to compute the full Mic scatter- 
ing properties of this dust model, the bhmie routine from 
(Bohren & Huffman 1983) and modified by B. Draine® was 
used. 

The radiative equilibrium temperatures were first com- 
puted, followed by images at 12 wavelengths ranging from 
1.25 microns to 500 microns^. The wavelengths were cho- 
sen to produce synthetic images in common bands, in- 
cluding JHK (1.25, 1.65, and 2.15 /im), Spitzer /IRAC 3.6, 
4.5, and 8.0 //m, Herschel /PACS (70, 110, and 170 ^m), 
and Herschd/SFIRE (250, 350, and 500 /im). A distance 
of 300 pc was used to simulate a nearby low-mass star- 
formation region. The images were computed at a viewing 
angle of {6, (/>) = (45°, 45°), and are shown in Figure 8 at a 
common resolution of l"/pix. 

The near-infrared (JHK) image shows emission from the 
diffuse cloud, which is scattered light from the interstel- 
lar radiation field, as well as bright features corresponding 
to light from the forming stars scattered in the circum- 
stcUar dust. The near-infrared images are entirely domi- 
nated by direct and scattered stellar light, and there is no 

^ The stellar component provides photon packets that are 
scattered off the cloud, but the direct (un-scattered) emission 
from this component is not included in the images presented in 
Figures 8 and 9 since it would not appear as a diffuse component 
but as point sources. 
^ http : //www. astro .princet on. edu/~ draine/ scattering.html 
^ The calculation used 10** photons for each specific energy 
iteration, lO' photons per wavelength for the monochromatic 
peeling-off for scattering for both the source and dust emission, 
and lO" photons for the source and dust emission raytracing. 
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Fig. 10. Spectral Energy Distributions for the three sourees labeUed in Figure 8. Eaeh panel shows SEDs for 20 different 
viewing angles. The black lines show the minimum and maximum values as a function of wavelength, the light gray shows 
the range of values, and the dark gray shows the individual SEDs. 



contribution from thermal emission. The mean colors of 
the near-infrared emission are J-H:= 0.387 ± 0.280 and H- 
K= —0.214 ± 0.273 where the uncertainties indicate the 
scatter in the values rather than the error in the mean. 

The IRAC image also shows (in blue) the scattered light 
from some of the sourees, and a couple of the sources show 
a bipolar structure. This is noteworthy because the simula- 
tion does not include outflows. The bipolar structures seen 
here are instead due to beaming of the radiation by the 
three-dimensional structure of the accretion flows and the 
presence of protostellar disks. These cause the density to be 
lower along the rotation axis, and the radiation to prefer- 
entially escape in those directions. Thus, bipolar scattered 
light structures in the near- and mid-infrared are not nec- 
essarily evidence for actual outflows. 

The IRAC image also shows the bright PAH back- 
ground from the interstellar radiation field, which causes 
the cloud to appear in absorption as an infrared-dark cloud 
(IRDC). The column density of the filaments is typical of 
IRDCs: for instance, 10% of the area shown in Figure 8 
has S > 0.10 g/cm^ or Ay > 34.6 mag and 1% of the area 
has E > 0.27 g/cm^ or Ay > 90.7 mag. In contrast, the 
far-infrared images show the cloud in emission. The color 
gradients across the PACS image trace gradients in the dust 
temperature, with the regions immediately surrounding the 
protostars being warmer than the rest of the cloud. Finally, 
the SPIRE image traces most of the mass, but does not 
show very strong color gradients, as most of the emission is 
in the Rayleigh- Jeans part of the spectrum. 

5.3. Synthetic Observations 

Synthetic observations were produced by resampling the 
pixel resolution of the images to values typical of ground 
based observations for JHK (1"), and of Spitzer and 
Herschel observations for the other bands (1.2" for IRAC, 
3.2" for PACS, and 6" for SPIRE). The observations were 
convolved with the appropriate instrumental PSFs and 
Gaussian noise with typical values (listed in Table 1) was 
added to the images. The resulting synthetic observations 
arc shown in Figure 9. 

The near-infrared (JHK) image does not change, since 
the pixel resolution is the same as before, and the PSF 



full-width at half maximum (1") is the same as the pixel 
resolution. The IRAC image shows many of the protostars 
appearing as strong 8 ^m point sources. These were not 
seen in Figure 8 since only one pixel contained this bright 
emission for each source. The PACS and SPIRE images 
are severely degraded, with only the brightest diffuse emis- 
sion seen in the PACS image. Nevertheless, the PACS im- 
age shows very strong compact emission associated with 
the protostars. The SPIRE image does show bright point 
sources for a few of the protostars, but the remainder are 
confused with the cloud. 

5.4. Spectral Energy Distributions 

For three sources, indicated in Figure 8, SEDs were com- 
puted inside 20 apertures, with a central circular aperture 
5,000 AU in radius, and the remaining apertures consist- 
ing of annuli logarithmically spaced between 5,000 AU to 
20,000 AU. The SEDs were directly computed by binning 
photon packets into apertures in the radiative transfer code 
(rather than using the images). In reality, measuring SEDs 
from convolved and noisy images is going to affect the re- 
sulting SEDs, but since this is highly dependent on the 
convolution and noise parameters, we do not take this into 
account here as we are interested in the 'intrinsic' SEDs. 
The SEDs were computed for 20 viewing angles regularly 
distributed in three dimensions around the sources*. The 
innermost circular aperture was used to compute the SED, 
and the median of the SEDs in the 19 surrounding annu- 
lar apertures was used as an estimate of the wavelength- 
dependent background that was subtracted from the cen- 
tral aperture SED. 

The resulting SEDs are shown in Figure 10. The three 
sources were chosen as they were expected to produce dif- 
ferent SEDs based on the images shown in Figure 8, and 
indeed the three sources show SEDs with varying amounts 
of near- and mid-infrared emission. The dependence of the 
near- and mid-infrared emission on viewing angle is very 

* Truly regular non-random three-dimensional spacing of 
points on a sphere is only possible for 4, 6, 8, 12, and 20 points, 
and the position of the points corresponds to the vertices of 
regular polyhedra. In the current section, the regular spacing is 
achieved by using the 20 vertices of a dodecahedron. 
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Table 2. Properties of the sources in Figure 10. 



Source 




M 


M/M-, 






^disk 


itot 




(Mq) 


(Afo-yr-i) 




Lo 


Le 


Le 


Le 


1 


2.216 


6.012e-07 


2.713e-07 


24.20 


1.25 


2.51 


27.96 


2a 


0.432 


3.192e-06 


7.388e-06 


0.03 


2.81 


5.62 


8.46 


2b 


1.548 


4.089e-06 


2.642e-06 


5.32 


9.74 


19.48 


34.54 


3 


0.058 


4.771e-06 


8.163e-05 


0.00 


1.45 


2.91 


4.36 



large, and is a consequence of the three-dimensional struc- 
ture of the circumstellar material. The basic evolutionary 
properties of the sources are listed in Table 2. Source 2 
is actually composed of a pair of sink particles, so these 
are listed as 2a and 2b in the table. Sources 1, 2a, and 2b 
have similar masses, of the order or 0.5 — 2Mq. Source 3 
in contrast has a much lower mass, around 0.06 Mq. By 
looking at the accretion rate of the sources relative to the 
sink particle masses, one can see that the relative accretion 
rate increases from Source 1 to 3, which is consistent with 
the average behavior of the near- to far-infrared ratio in 
the SEDs. However, the evolutionary stage of the sources - 
both relative to each other and in absolute terms - may be 
derived incorrectly in reality, when only one viewing angle 
is available. 

6. Future 

The code will be actively developed in the future to im- 
prove existing features and add new capabilities. Some of 
the planned capabilities include: 

— Scattering and absorption by non-spherical grains 
aligned along magnetic fields. The main effect of 
aligned grains is to produce linear and polarization 
polarization signatures that trace the magnetic fields 
(Whitney & Wolff 2002). 

— Continuum and line gas radiative transfer, including 
photoionization. This will be very useful for modeling 
ionization and dynamics in a wide range of environ- 
ments. 

— Support for unstructured meshes, where each cell is 
an irregular polyhedron, and which are commonly con- 
structed from a Voronoi tessellation of space. They allow 
a more precise sampling of arbitrary three-dimensional 
density structures than other types of adaptive grids 
for a fixed number of cells. Unstructured meshes are 
now being used in hydrodynamical simulations (e.g. 
Springel 2011; Greif et al. 2011), and support for these 
in Hyperion would allow radiative transfer post- 
processing of the simulations without the need for re- 
gridding the density structures. 

— Ray tracing for scattered light, which - while memory 
intensive - would provide high signal-to-noise for both 
images and SEDs efficiently at wavelengths dominated 
by scattered light faster than the peeling-ofF algorithm 
(Pinte et al. 2009). 

— Temperature dependent opacities for dust, since the 
composition of dust often depends on the temperature, 
especially when including ices. 

7. Summary 

Hyperion is a new radiative transfer code that has been 
designed to be applicable to a wide range of problems. 



Radiative transfer can be carried out through a variety 
of three-dimensional grids, including cartesian, cylindrical 
and spherical polar, and adaptive cartesian grids. The al- 
gorithms used by the code, including the photon packet 
propagation, raytracing, convergence detection, and diffu- 
sion approximations are described in detail (§2). 

The code is parallelized and scales well to hundreds or 
thousands of processes. As expected, the maximum speedup 
depends on the fraction of the execution time spent in the 
serial part of the code (such as file input /output), which in 
turn depends on the problem being computed (§3). 

The code was tested against the Pascucci et al. (2004) 
and Pinte et al. (2009) benchmarks for circumstellar disks 
and was found to be in excellent agreement with the re- 
spective published benchmark results (§4). 

As a case study, Hyperion was used to predict near- 
infrared to sub-millimeter observations of a star formation 
region using a dynamical simulation of a region of low-mass 
star formation (Offner et al. 2009). One interesting finding 
was that despite the absence of outflows in the simulation, 
three-dimensional accretion flows caused beaming of radia- 
tion that produced fairly symmetric structures usually as- 
sociated with bipolar outflow cavities. 

Hyperion is published under an open-source license 
at http://www.hyperion-rt.org, using a version control 
hosting website that makes it easy for users to contribute 
to the code and documentation. 
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Appendix A: The partial difFusion approximation (PDA) 

The diffusion approximation is described by 

V • (DVT'^) = 0, 



(A.l) 



where D — l/3pxR is the diffusion coefficient. This equation can be written as a system of hnear equations. Since 
Hyperion stores the specific energy absorption rate of the dust rather than the temperature, the PDA is actuaUy solved 
using: 



DV 



A 



Kp{A) 



0. 



In spherical polar coordinates, the diffusion equation is 



j,2 Qj. 



r^D 



dr 



sin 6 rdO \ rdO 



d 



r sin 9d(j) 



D 



r sm 



= 0. 



This equation can be written for a discrete spherical polar grid as 
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= 0, 



(A.2) 



(A.3) 



(A.4) 



where i, j, and k are the indices of the grid cells, r, and are the coordinates of the center of the cell, Ar^, rA9j, and 
7'sin6'A(/)fc are the spatial extent of the cell along each coordinate, and the Ar terms are the Rosseland optical depth 
across the cell along each coordinate. This equation can be written down for all the cells in which the PDA is required. 
In some cases, the T value will be well determined if they neighbor the diffusion region - these values act as boundary 
conditions. The result is a system of linear equations that can be formally solved. 
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